NASA Technical Memorandum 106258 


/// - 

m 


Desktop Chaotic Systems: Intuition 
and Visualization _ 


Michelle M, Bright and Kevin J. Melcher 

Lewis Research Center 

Cleveland, Ohio 


and 


Helen K. Qammar and Tom J'. Hartley 

The University of Akron 

Akron, Ohio z " 


(NASA-TM-106253) OESKTOP CHAOTIC 
SYSTEMS: INTUITION AND 

VISUALIZATION (NASA) 14 p 


N93-31847 


Unci as 


G3/64 0175559 


Prepared for the" 

SPIE International Symposium on Chaos/Nonlinear Dynamics 
sponsored by Society of Photo-Optical Instrumentation Engineers 
San Diego, California, July 11-16, 1993 





Desktop Chaotic Systems: Intuition and Visualization 

Michelle M. Bright and Kevin J. Melcher 
NASA Lewis Research Center 

Advanced Controls Technology and System Dynamics Branches 
Cleveland, OH 44135 

Helen K. Qammar and Tom T. Hartley 

Departments of Chemical and Electrical Engineering 
The University of Akron 
Akron, Ohio 44325 


Abstract 

This paper presents a dynamic study of the Wildwood 0 Pendulum, a commercially available 
desktop system which exhibits a strange attractor. The purpose of studying this chaotic pendulum is two- 
fold: to gain insight in the paradigmatic approach of modeling, simulating, and determining chaos in 
nonlinear systems, and to provide a desktop model of chaos as a visual tool. For this study the nonlinear 
behavior of this chaotic pendulum is modeled, a computer simulation is performed, * and an experimental 
performance is measured. An assessment of the pendulum in the phase plane shows the strange attractor. 
Through the use of a box-assisted correlation dimension methodology, the attractor dimension is 
determined for both the model and the experimental pendulum systems. Correlation dimension results 
indicate that the pendulum and the model are chaotic and their fractal dimensions are similar. 

2. Motivation 

When studying nonlinear dynamics and chaotic behavior, there is a need to clearly present the 
concepts of chaos, strange attractors, and correlation dimension to a non specialist. Many software tools 
exist to assist in the graphical display and visualization of chaotic phenomena. However, there is a 
definite need for a variety of desktop chaotic systems, that is, small scale devices that clearly demonstrate 
chaotic behavior to a variety of audiences, yet are easily modeled, and can be used as teaching tools. 

Some of these desktop chaotic devices have been identified in Moon [1], for example, the 
coupled neon bulb chaotic relaxation oscillator and Rott’s coupled pendula. "The neon bulb system 
consists of two neon bulb circuits coupled together and the two circuits exhibit chaotic dynamics in the 
form of the flashing bulbs." The other example, Rott’s coupled pendula, "allows for transient chaotic 
behavior and is designed like a three-armed pendulum puppet". Both these systems are widely exercised 
to gain intuition into their dynamics and for a means to visualize chaos. Other constructable desktop 
systems can be found in [1]. In this paper another chaotic desktop system is presented, the Wildwood 
Pendulum 0 This system is chosen because it exhibits chaotic motion, is commercially available, is 
inexpensive, and its unpredictable behavior is intuitive and easily observed. 
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3. System Description 


The Wildwood Pendulum 0 is a ferromagnetic end mass on a flexible shaft oscillating above an 
electromagnet and five permanent magnets. The pivot point is a magnetic hinge with very low damping. 
Figure la shows an isometric drawing of the pendulum. Figures lb and 1c provide top and side views 
which show the arrangement of the pendulum system and define key parameters. Table A provides a list 
of measured parameter values used in the pendulum model. 


PARAMETER 

DESCRIPTION 

VALUE 

UNITS 

G 

gravitational constant 

9.8 

m/s 

w 

width of base 

93 

mm 

1 

length of base 

164 

mm 

M 

mass of pendulum (rod and bob) 

6.935 

grams 

C 

pendulum damping factor 

0.015 


k P 

field strength of bob magnet 

300 

gauss 

R, 

length of pendulum bob 

145 

mm 

R„ 

distance from hinge to base plane 

150 

mm 

N ra 

number of permanent magnets in base 

5 


K 

radius of permanent magnets in base 

26.77 

mm 

r m 

radius of each permanent magnet 

6.35 

mm 

k m 

field strength of permanent magnet 

100 

gauss 

Re 

radius of electromagnet 

8.5 

mm 


Table A - Parameter Values for the Wildwood Pendulum 0 


The pendulum’s chaotic behavior is observed when the bob is perturbed about the unforced 

equilibrium position, Q o , which is shown in Figure lc. As the bob moves toward 0 o , a sense coil 

engages the electromagnet which exerts an attractive force on the bob. As the bob moves away from 

0^ , the sense coil disengages the electromagnet. This has the effect of adding energy to the system by 

increasing the acceleration of the pendulum bob. A diagram of the patented circuit for the electromagnet 
and sense coil is shown in Figure 2. While the electromagnet exerts an attractive force on the pendulum, 
the ferromagnets exert a constant repulsive force. 

4. Equations of Motion 

To better understand the dynamics of this pendulum and to improve our knowledge of chaotic 
systems, a mathematical analysis of this system is performed. To model the pendulum, two different 
coordinate systems are used for ease of implementation, and to relate the magnetic forces of the 
pendulum bob to the magnetic forces located in the pendulum base. The coordinate system for the base 
is shown in Figure 3a. The base plane coordinates have unit vectors e x , e y , and e z , in the directions of the 
x, y, and z axes respectively. Here, the origin is located in a horizontal plane parallel to the top of the 
permanent magnets and directly beneath the pendulum hinge point. Spherical coordinates for the 
pendulum bob are shown in Figure 3b. This coordinate system has unit vectors e e , e^, and e r in the theta, 
phi and radial directions. The origin of the spherical coordinate system is at the pendulum hinge 

point, 0 o . 
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In order to simplify the model, the pendulum shaft is assumed to be rigid and of fixed length. 
Therefore, r is a constant and velocity, f , and acceleration, r , in the radial direction are zero. In this 
case the spherical equations of motion [2] reduce to: 

position: r = re r 


' velocity: 
acceleration: 


v = r0e e + r4>sin0e. 


a= (r0 -r^sinGcosG)^ 


0 * 


m 


2 4> A 
m 


(r$sin0 +2r0<j)cos0)e <{) 


Solving the acceleration equation for the acceleration terms 0 and <$>, yields the following 
differential equations which may be integrated to obtain the pendulum’s velocity and position. 


V 


— + r<j) 2 sm0cos0 
m 


V 4> 


1 

A 

rsinG 

m 


- 2r0<j>cos0 


To solve for the required derivatives, the forces exerted on the pendulum bob by the various 
magnets must be computed. For simplicity, it is assumed that the force of a particular magnet on the bob 
is proportional to the inverse of the distance squared between the pendulum and the magnet. The force 
equation is then: 

KK KK 

Y — p m — p m 

2 2 2 2 
Ar 2 AX +*y +AZ 

Here, K p is the force constant for the permanent magnet in the pendulum bob, and K m is the force 
constant for permanent magnets in the base. The distance between the two magnets in question, at, is 
computed using base plane coordinates. The sign of the magnetic force constants determines whether the 
force is attractive or repulsive in nature. If signs are opposite, the magnets attract. If signs are the same, 
they repel. 

This same method is used to compute the force of the electromagnet on the pendulum bob. As 
the magnetic pendulum bob accelerates towards the electromagnet, the above force equation is assumed; 
however, when the pendulum bob moves away from the electromagnet, the force is assumed to be zero 

because the electromagnet is switched off when 0 < 0. 

When all of the forces have been computed, forces in the same coordinate direction are summed, 
the resulting base plane force vector is then transformed to the spherical coordinate system of the 

pendulum bob. These forces are then used to compute the acceleration terms 0 and 
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At this point, some discussion concerning the peculiarities of the model is in order. Ail of the 
parameters were measured, except for the magnetic strength of the electromagnet. This parameter is 
given a relative magnitude with respect to the measured values of the other magnets and observed 
behavior of the pendulum. Initially, an attempt was made to measure this value, however, due to the 
gaussmeter used, and the interacting field strengths of the base magnets, an accurate measurement of the 
electromagnet was not made. 

Additionally, a single point source was used initially to represent each magnet. This was 
insufficient due to the interacting field strengths of the magnets, because the size of the base magnets is 
the same order of magnitude as the distance between the magnets. Representing each magnet in the base 
as a cluster of point sources produced behavior more accurately representative of the physical system. 

This clustering of point sources is shown in Figure 4. 

The simulation of the Wildwood Pendulum 0 was performed using the MATLAB^ 4.0 matrix 
math package. A copy of the simulation code is given at the end of this paper. Figure 4 gives a typical 
response of the simulation in the x-y plane. This time response of the system was obtained by integrating 
values of the differential equations using a Runge-Kutta-Fehlberg method. Figure 5 shows a phase plane 

representation of this simulation data for 0 vs. 0. This phase plane representation gives evidence of the 
chaotic nature of this pendulum simulation. This simulation data is used in the dimension determination 
section to calculate an estimate of the correlation dimension. 

5. Experimental Data 

To analyze the chaotic behavior of the actual pendulum system, data was gathered from the 
Wildwood Pendulum 0 using the photodiode circuit shown in Figure 6. As the pendulum passes through 
the beam of light, the voltage across the photodiode decreases in proportion to the amount of light 
observed. This is equivalent to capturing the dynamics of the pendulum in the plane of the light beam. 
The time series plot of the voltage over time for the pendulum is shown in Figure 7. The time between 
zero crossings of the voltage signal represent the chaotic nature of the pendulum system. 

To compare this voltage data with the simulation data, a phase plane portrait of the voltage data 
is reconstructed using 10000 data points. To reconstruct the attractor, an independent time series of data 
was generated using the following time delay embedding method[3]: 

v(t) = (x(t), x(t+T), ..., x(t+(m-l)T). 

where v(t) is the independent time series, x(t) is the voltage signal, m is the embedding dimension, and T 
is the reconstruction time delay. For this voltage data, the values of m=4 and T=1 were used for the 
reconstruction. This reconstruction is used in the next section for determining the dimension of the 
experimental system. 


6. Dimension Determination 

Verification of the fractal behavior of this pendulum is provided using a box-assisted correlation 
dimension algorithm [4,5] to determine the dimension of both the simulation data and the measured 
pendulum data. Two box-counting correlation dimension algorithms are used to ensure the dimension is 
accurate. 
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For the correlation dimension method, 4000 data points are used and the first 1000 points are 
thrown out as transient. Using 3000 data points for both the experimental data and the simulation data 
allows for accurate prediction of correlation dimension up to approximately 7, according to the following 
equation from Eckmann and Ruelle[6] which states: 

Maximum Accurate Dimension < 2 log(number of data points). 

With this limit in mind, the correlation dimension, (or correlation exponent), was calculated for both 
pendulum systems using a Grassberger-Procaccia correlation integral. [7] This correlation exponent is 
closely related to the fractal dimension and gives a lower bound to the fractal dimension. The 
Grassberger-Procaccia algorithm using Theiler’s method[4,5] considers the spatial correlation between 
pairs of points on a reconstructed attractor, and is measured with the following correlation integral: 

C(N,r) = _2_ I H( r - || Xi - Xj ||), 

N(N - 1) 

where H(x) is the Heaviside step function. The summation counts the number of pairs ( x„ Xj ) for which 
this distance is less than r, where r is the box size. Then C(N,r) scales like a power of v so that: 

C(N,r) = r v 

where v is the correlation dimension which is the slope of the log-log plot of C(N,r) vs. r. 

The correlation dimension is calculated first for the experimental pendulum data with Theiler’s 
box-assisted method, using embedding dimensions from m=4 to m=7. The correlation dimension from 
this method is approximately 3.47. The same correlation dimension method is used to calculate the 
dimension for the simulation data. The log-log plots of C(N,r) vs. r for the simulation data only 
converged for embedding dimensions of m=4 and m=5, with a resulting* correlation dimension of 3.2. 

To corroborate the dimension results for the simulation data, another box-counting method using 
the Grassberger-Procaccia algorithm was implemented. This method was used only on the simulation 
data. From this analysis the simulation equations exhibit a fractal dimension of approximately 3.3. 

From these dimension results it is concluded that the modeled pendulum equations closely match 
the behavior exhibited by the actual Wildwood 0 chaotic pendulum. Both the simulated data and the 
measured data resulted in a similar correlation dimension of approximately 3.4. 

7. Conclusions 

From the study of desktop chaotic systems insight is gained in modeling, simulating, and 
analyzing chaos in nonlinear systems. The Wildwood ° pendulum is chosen for analysis because its 
chaotic behavior is easily observed and its fractal dimension has not been determined previously. For this 
pendulum, the model equations are given and analyzed using both correlation dimension techniques and 
Lyapunov characteristic exponent methods. Experimental data is gathered for the pendulum, and the 
attractor dimension of the measured data is calculated. Both the modeled pendulum and the measured 
pendulum data represent a chaotic system with similar attractor dimension. This process was made easier 
with the use of desktop systems like the chaotic pendulum, which provided a system to easily visualize 
chaotic behavior. 
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Figure 1c. Side View of the Chaotic Pendulum Figure 2. The Electromagnet and Sense Coil Circuit 



Figure 3a. Base Plane Coordinate System Figure 3b. Spherical Coordinate System 
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Theta'(t), radlans/sec 




light source 


photodiode 


Figure 4. X-Y Plane Plot of Pendulum Simulation 


Figure 6. Photodiode Circuit for Experiment 


CHAOTIC PENDULUM SIMULATION 



Theta(t), radians 


Figure 5. Phase Plane Plot of the Pendulum Simulation Figure 7. Voltage Across Photodiode Versus Time 
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%PEND3RUN.M ... Matlab M-file used to Run Wildwood Pendulum Simulation 


%Global Variables ... 
global G\ 

% global W; 
global M; 
global Kd; 
global Kp; 
global Rp; 
global Rb; 


^Gravitational Cosntant 
^Weight of Pendulum 
%Mass of Pendulum 
%Pendulum Damping Constant 

%Pendulum Magnetic field Strength Const (Perm Magnet) 
%Pendulum Radius 

^Distance from Pendulum Support to Base Plane 


% global Nm; 

% global Mm; 
% global NMm; 
% global Rm; 

% global rm; 
global Xm; 
global Yra; 
global Km; 


%Number of Permanent Magnets 
%Number of Point Sources/Perm Magnet 
%Total Number of Point Sources for Perm Magnets 
%Radius of Permanent Magnets in Base Plane 
%Radius of Individual Permanent Magnet 

%X Loc of Permanet Magnet Sources in Base Plane Coords 
%Y Loc of Permanet Magnet Sources in Base Plane Coords 
%Field Strength Const of Individual Perm Magnets) 


% global NMe; 
% global re; 
global Xe; 
global Ye; 
global Ke; 


%Total Number of Point Sources for Electro Magnet 
%Radius of the Electro Magnet 

%X Loc of Electromagnet Sources in Base Plane Coords 
%Y Loc of Electromagnet Sources in Base Plane Coords 
%Field Strength Const of Electromagnet 


%If "CONTINUE" is set then start simulation from where it left off ... 
if exist(’CONiJLNUE’) ~ = 1; 

%K0 is used to Vary different parameters in the simulation 
ifexist(’KO’)~==l; K0=1, end; 


%Define Parameters for Pendulum Bob 
G = 9.8; %meters/sec/sec 

M = 0.007; %kg 

W = M*G; %newtons 

Kp = 3.0e-3; 

Rp = 0.145; %meters 

Rb = 0.150; %meters 


%C alculate some reasonable numbers for damping 
% Damping = -Kd*THd 

zeta = 0.015; %daraping constant for pendulum 

wn = sqrt(G/Rp); %undamped natural frequency of simple pendulum 

Kd = 2*zeta*wn; %kg/sec 

%Define Parameters for Permanent Magnets 
Rm = 0.027; %meters 

rm = 0.007; %meters 

%Use 5 Evenly Spaced Magnets. Each Magnet has 7 Point Sources. 

Nm = 5; %Number of Permanent Magnets (best odd) 

Mm = 7; %Nuraber of Point Sources/Magnet (even perimeter best) 

NMm = Nm*Mm; %TotaI Number of Magnets 

^Calculate Position of Magnets in Base Plane Coords 
THc = 2*pi*[0:Nm-l]/Nm; %Theta Position Each Magnet 

%Convert Magnet Position From Polar to Cart Coords 
[Xc,Yc] = pol 2cart(THc,Rra * ones(size(THc))); 

^Calculate Position of Point Sources relative to center of each magnet 
% Save one point source for the center 
THm = 2*pi*[0:Mm-2]7(Mm-l)+pi; %Theta Position 

%Convert Point Source Position to Cart Coords 
[Xt, Y t] = po!2cart(THra ,rm * ones(size(THm))); 

%Source Absolute Position = Source Relative Position + Magenet Mid Point 
% Set Source at Mid Point of each magnet 
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Xni = [Xc; ones(sizc(Xt))*Xc+Xt*cos(THc)-Yt*sin(THc)]; Xra=Xra(:y; 
Yra = [Yc; ones(size(Tt))*Yc+Xt*sin(THc)+Yt*cos(THc)]; Yra=Ym(:y; 
%Forcc Const for Each Magnet 
% Need same sign as Kp for repulsive force 
Km = (Kp/3)*ones(size(Xm))/Mra; 

%Define Parameters for Electro Magnets 
re = 0.008; %meters 

%Use a Single Magnet with 9 Point Sources. 

NMe = 9; %Total Number of Magnets 

^Calculate Position of Point Sources relative to center of magnet 
% Save one Point Source for Center of each magnet 
THe=2*pi*[0:NMe-2] 7(NMe-l)+pi; %Theta Position 

^Convert Point Source Position to Cart Coords 
pCe, Ye] = pol2cart(THe t re * ones(size (THe))); 

%Set Point Source at center of each magnet 

Xe=[0 Xe(:y]; 

Ye=[0 Ye(:)’]; 

%Force Const for the Magnet 
% Need sign opposite of Kp for attraction 
Ke = -Kp/40*ones(size(Xe))/NMe; 

%Define Start (T) and Stop (Tfinal) Times for Integration Routine 
T=0; Tfinal = 10; Tc2Ts=0; 

9&Define Initial Conditions of the Pendulum Relative to Hinge Point 
%... [Theta ThetaDOT Phi PhiDOT] (angles in radians) 

% 0 >= Theta <= pi 

% -pi >= Phi <= pi 
STATES = [pi-pi/ll 0 pi/4 0]; 
end; 

%Break out of Integration Routine to save data and report to user 
% every DTR seconds 
DTR= (Tfinal -m ax (T))/20; 

%Run System Simulation 
pend2plt; set(gca, , DrawMode , /fast’); drawnow 
if ( max(T)<Tfinal-eps ) 

TcO = clock; TsO=raax(T); 
while (max(T)+eps<Tfinal) 

^Calculate Time Limits for Current Segment of Simulation 
n=Iength(T); tO=T(n); tl=min(tO-i-DTRTCnal); 

^Report Status to the User 

note = [’Running the Simulation from T=’,n2s(t0),’ to T=’,n2s(tl)]; 
if Tc2Ts>0; note=[note Qock/Sim TIrae=’,n2s(Tc2Ts)]; end 
disp(note) 

%Run Simulation for Specified Segment 
[t,states] = ode23(’pend3ode\t0,tl,STATES(n,:) t 2.0e-4); 

% Append Results to Existing States 

[n,m]=size(states); T=[T;t(2:n)]; STATES=[STATES;states(2:n,:)]; 
%PIot Current Data 
hold on; 

pIot(Rp*sin(states(:,l)).*cos(states(:,3)), ... 
Rp*sin(states(:,l)).*sin(states(: T 3)),’rn.\ ... 

TraseMode ’/n on e ’ ) ; 
hold off; 

if exist(’hTf)= = l; 

set(hTf,’String’, [Tfinal = ’^um2str(max(T)),’sec’]); 
end; 

drawnow 

%Give User a Measure of the System Performance 
Tc2Ts = fix(etime(clock T TcO)/(max(t)-TsO)); 
end 
end; 
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function [derivs] = pend2odc(t, states) 

% 

% FUNCTION: pend2ode 

% 

% PURPOSE: Compute derivatives for a conical pendulum in a chaotic system. 
% 

% SYNTAX: [derivs] =pend2ode(t, states) 

% 


% Global Variables ... 

global G; 

% global W; 
global M; 
global Kd; 
global Kp; 
global Rp; 
global Rb; 


% Gravitational Cosntant 
%Weight of Pendulum 
%Mass of Pendulum 
%Pendulum Damping Constant 

%Pendulum Magnetic Geld Strength Const (Perm Magnet) 
%Pendulura Radius 

%Distance from Pendulum Support to Base Plane 


% global Nm; 

% global Mm; 
% global NMm; 
% global Rm; 

% global rm; 
global Xm; 
global Yro; 
global Km; 


%Number of Permanent Magnets 
%Number of Point Sources/Perra Magnet 
%Total Number of Point Sources for Perm Magnets 
%Radius of Permanent Magnets in Base Plane 
%Radius of Individual Permanent Magnet 

%X Location of Permanet Magnet Sources in Base Plane Coords 
%Y Location of Permanet Magnet Sources in Base Plane Coords 
%Field Strength Const of Individual Perm Magnets) 


% global NMe; 
% global re; 
global Xe; 
global Ye; 
global Ke; 


%Total Number of Point Sources for Electro Magnet 
%Radius of the Electro Magnet 

%X Location of Electromagnet Sources in Base Plane Coords 
%Y Location of Electromagnet Sources in Base Plane Coords 
%Field Strength Const of Electromagnet 


%Extract States from "state” vector 

%THETA: measured counter- cl ockwise from the x-axis in the x-y plane) 
TH = states(l); %Theta 
THd = states(2); %Theta dot 

%PHI: measured clockwise from the z-axis in the y-z plane) 

PH = states(3); %Phi 
PHd = states(4); %Phi dot 

%Calculate Location of Pendulum in Cart Base Plane Coords 
Xp = Rp*sin(TH)*cos(PH); 

Yp = Rp*sin(TH)*sin(PH); 

Zp = Rb+Rp*cos(TH); 


^Calculate the Distance Between the 
dXm = Xp*ones(size(Xm))-Xm;% 
dYm = Yp*ones(size(Ym))-Ym;% 
rpsq = dXm.^2 + dYm. ^ 2;% 
rp = sqrt(rpsq);% 
rmsq — rpsq + Zp*Zp;% 
rm = sqrt(rmsq);% 


Permanent Magnets and the Pendulum 

%x dist pend to mag 

9cy dist pend to mag 

%dist in plane squared 

%dist in plane 

%dist pt to pt squared 

%dist pt to pt 


^Calculate Forces of Permanent Magnets on Pendulum 
% Assume that the Force is proportional to Kl*K2/(r*r), 

% where Y is the distance between the pendulum and a magnet 
% If Kp*Km < 0 force is attracting, 

% If Kp*Km = 0 No force 
% If Kp*Km > 0 force repelling 
Fm = Kp*Km./rmsq;% %Total force 
%Reduce the Forces to Their Cart Coord Components 
Fmx = Fm.*dXm./rm;% %x force in base plane coords 

Fmy = Fm.*dYm./rm;% %y force in base plane coords 
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Fraz = Fm.*Zp./rm;% 


%z force in base plane coords 


96Ca3uculate Force of Electromagnet cluster (at 0;0) on pendulum 
% Assume that the Force is proportional to Kp*Ke/(r*r), 

% where V is the distance between the pendulum and a electromagnet 
% If Kp*Ke < 0 force is attracting, 

% If Kp*Ke = 0 No force t 

% If Kp*Ke > 0 force repelling 

dXe = Xp*ones(size(Xe))-Xe; %x dist pend to mag 

dYe = Yp*ones(size(Ye))-Ye; %y dist pend to mag 

rpsq = dXe. ~2 + dYe. ^2; %dist in plane squared 

rp » sqrt(rpsq); %dist in plane 

resq = rpsq + Zp*Zp; %dist pt to pt squared 

re = sqrt(resq); %dist pt to pt 

%The Electro Magnet is on only if the Pendulum Bob is moving 
% toward it 
if ( THd >0) 

Fe = Kp.*Ke./resq; %Total force 

else 

Fe = 0; 
end 

^Reduce the Forces to Their Cart Coord Components 
Fex = Fe.*dXe./re;% %x force in base plane coords 

Fey as Fe.*dYe./re;% %y force in base plane coords 

Fez = Fe.*Zp./re;% %z force in base plane coords 

^Calculate Gravitational Force on Pendulum 
Fgz = -M *Gr y % %Force Due to Gravity in base plane coords 

^Calculate Total Forces in Base Plane 

Ftx = sum( [ Fmx, Fex ] )\% %x Force in base plane coords 

Fty = sum( [ Fmy, Fey ] );% %y Force in base plane coords 

Ftz - sum( [ Fmz, Fez, Fgz ] );% %z Force in base plane coords 

^Transform Base Plane Forces to Spherical Coords 

% Define Transformation Matrix 
TM=[ cos(PH) * cos(TH) sin(PH)*cos(TH) -sin(TH) 

-sin(PH) cos(PH) 0 
cos(PH)*sin(TH) rin(PH)*sin(TH) cos(TH) 

]; 

^Calculate the Forces in the Pendulums Coord 
FTH = TM(1 , :) * [Ftx;Fty;Ftz] ; 

FPH = TM(2,:)* [Ftx;Fty;Ftz] ; 

FR = TM(3,:)*[Ftx;Fty;Ftz]; 

^Calculate the Derivatives 
RpSq = Rp*Rp; THdc = THd ^ 3; 

THdd = (Fin/M + Rp*PHd*PHd*sin(TH)*cos(TH) -Kd * THd * Rp)/Rp; 
if ( TH ~ = 0 ) 

PHdd = (FPH/M - 2*Rp*THd*PHd*cos(TH))/(Rp*sin(TH)); 
end 

derivs = [ THd THdd PHd PHdd]; 
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